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FRACTALS  AND  PHYSIOLOGY 


■  Fractal  characteristics: 


■  Complex,  patterned 

■  Statistical  self-similarity 

■  Scale-invariant  structure 

■  Generated  by  simple 
iterative  rules 

■  Mu>2H+d  spectral  decay 


QuickTime™  and  a 
TIFF  (Uncompressed)  decompressor 
are  needed  to  see  this  picture. 


■  Growth  processes,  biofractals 


QuickTime™  and  a 
TIFF  (Uncompressed)  decompressor 
are  needed  to  see  this  picture. 


QuickTime™  and  a 
TIFF  (LZW)  decompressor 
are  needed  to  see  this  picture. 


QuickTime™  and  a 
TIFF  (Uncompressed)  decompressor 
are  needed  to  see  this  picture. 


QuickTime™  and  a 
TIFF  (Uncompressed)  decompressor 
are  needed  to  see  this  picture. 


Cardiovascular  system 


■  Heart 

■  Arterial  tree 

■  Dendritic  anatomy 


■  Lung 


QuickTime™  and  a 
TIFF  (LZW)  decompressor 
are  needed  to  see  this  picture. 


QuickTime™  and  a 
TIFF  (LZW)  decompressor 
are  needed  to  see  this  picture. 


QuickTime™  and  a 
TIFF  (LZW)  decompressor 
are  needed  to  see  this  picture. 


From  Goldberger,  Rigney  and  West 
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Fractal  bones 


■  Trabecular  bone 


QuickTime™  and  a 
TIFF  (LZW)  decompressor 
are  needed  to  see  this  picture. 


CT  of  a  vertebra 


__  QuickTime™  and  a 
TIFF  (LZW)  decompressor 
are  needed  to  see  this  picture. 


(jCT 


Courtesy  F.  Peyrin  ESRF 
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Mammograms 


QuickTime™  and  a 
TIFF  (Uncompressed)  decompressor 
are  needed  to  see  this  picture. 


QuickTime™  and  a 
TIFF  (Uncompressed)  decompressor 
are  needed  to  see  this  picture. 


QuickTime™  and  a 
TIFF  (Uncompressed)  decompressor 
are  needed  to  see  this  picture. 


QuickTime™  and  a 
TIFF  (Uncompressed)  decompressor 
are  needed  to  see  this  picture. 


DDSM:  University  of  Florida 

(Digital  Database  for  Screening  Mammography) 


(Arneodo  et  al.,  2001) 
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Brain  as  a  biofractal 


(Bullmore,  1994) 


QuickTime™  and  a 
TIFF  (LZW)  decompressor 
are  needed  to  see  this  picture. 


1mm 


Courtesy  R.  Mueller  ETHZ 
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OUTLINE 


■  Fractals  in  physiology  fp 

■  Wavelets  and  fractals 

■  Motivation  for  using  wavelets 

■  Fractal  processing:  order  is  the  key 

■  What  about  fractional  differentiation 

■  Fractional  splines 

■  Fractional  wavelets 

■  Wavelets  in  medical  imaging 

■  Survey  of  applications 

■  Analysis  of  functional  images 


Motivation  for  using  wavelets 


■  Wavelets  provide  basis  functions 
that  are  self-similar  [Mallat,  1989] 


V/O)  e  L2,  f(x)  =  ££</,&,*>  VSi,k(x) 


•  Wavelets  approximately  decorrelate  statistically  self-similar 
processes  [Flandrin,  1992;  Wornell,  1993] 

■  Unlike  Fourier  exponentials,  wavelets  are  jointly  localized  in 
space  and  frequency 

■  The  basis  functions  themselves  are  fractals  [Blu-Unser,  2002] 


Wavelets  are  prime  candidates  for  processing 
fractal-like  signals  and  images 
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On  the  fractal  nature  of  wavelets 


■  Harmonic  spline  decomposition  of  wavelets 

Theorem:  Any  valid  compactly  supported  scaling  function  (p{ x)  (or  wavelet  yA,x)) 
can  be  expressed  either  as 

(1)  a  weighed  sum  of  the  integer  shifts  of  a  self-similar  function  (fractal) ; 

(2)  a  linear  combination  of  harmonic  splines  with  complex  exponents. 

[Blu-Unser,  2002] 


QuickTime™  and  a 
TIFF  (Uncompressed)  decompressor 
are  needed  to  see  this  picture. 
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D4  as  a  sum  of  harmonic  splines 


■  Sum  of  spline  components 

+N/2 

<PN(X)  =  X  rnSn(x) 

n=-N/2 


where 


log  A  ,2nn 

*„(*)=  X/4(*-£)!082  7 1082 


k&Z  + 


QuickTime™  and  a 
Video  decompressor 
are  needed  to  see  this  picture. 
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Fractal  processing:  order  is  the  key  ! 


■  Vanishing  moments 

Classical  A/th  order  transform  <»  analysis  wavelet  y/{x)  has  N  vanishing  moments 


J  xnyA.x)dx  =  0,  n  =  0,---,N-l 

xs  R 


y/  kills  all  polynomials  of  degree  n  <N 


Multi-scale  differentiation  property 


Property 

An  analysis  wavelet  of  order  N 
acts  like  a  A/th  order  differentiator: 

yA^co)  =  0(coN ) 


□)  (fix),  y/(x-u))  =  -^{</>*f}(u) 

Smoothing  kernel:  (j){co )  =  y/  ( cl>)/(Jco ) 


CD 
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What  about  fractional  differentiation  ? 


■  Fractional  differentiation  operator 

^  fix)  <  F  >  {jcoj  f{a>)  yeR~ 


■  Motivation:  whitening  of  fBM-like  processes 

(j)A(o)  ~  0(1  / o)2H+l )  *  0(1) 

Fractional 

differentiator 

White  noise 

y=H+\ 


.  QUESTION 

Are  there  wavelets  that  act  like  fractional  differentiators  ? 

.  ANSWER 

Not  within  the  context  of  standard  wavelet  theory  where  the 
order  is  constrained  to  be  an  integer, 
but ... 
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SPLINES 


■  Polynomial  splines 

■  Fractional  B-splines 

■  Properties 

■  Fractional  differentiation 

■  Fractional  order  of  approximation 
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Polynomial  splines  (Schoenberg,  1946) 


Definition: 

s(x)  is  cardinal  polynomial  spline  of  degree  n  iff 

■  Piecewise  polynomial: 

s(x)  is  a  polynomial  of  degree  n  in  each  interval  [k,k+ 1); 

■  Higher-order  continuity: 

s(x),s{1\x),. .  .,s{n~l) (x)  are  continuous  at  the  knots  k. 
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B-spline  representation 


■  B-splines  of  degree  n 

(n+ 1)  times 

Explicit  formula  :  /3n+(x)  =  —  — 

n\ 


Theorem  [Schoenberg,  1946] 

A  cardinal  spline  of  degree  n  has  a  stable,  unique  representation 
as  a  linear  combination  of  shifted  B-splines 

s(x)=YJc(k)j3,:(x-k) 

k&Z 


Can  we  fractionalize  splines  ? 


■  Schoenberg’s  formula 


/?(*)  = 


n\ 


Pa+(x)  = 


a 

+ 


T(a  +  l) 
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Basic  tools  for  fractionalization 


■  Generalized  factorials — Euler’s  Gamma  function 

n\=  T(n  + 1)  Hw)  =  J  xu~'e~xdx 


Generalized  binomial 

+00 

(l  +  z/=X 


\ 

y 

k 

fwl 

UJ 

Z 

Ivy 

T(u  + 1) 


r(v  +  i)r(M-v  +  i) 


Fractional  derivative  [Liouvillle,  1855] 

&  <  Faaier  >  (jcoy 


Fractional  finite  differences 

^  Fourier  ^  ^  _e~j®y 


AVW  =  £(-!)' 


k=0 


(  > 
s 


V'V 


fix  -  k) 
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Fractional  B-splines 


/oO/x  0  /  ixO  y  Fourier 

J3+(x):=x+-(x-  1)+  < - > 


f  i  -jco  \ 

1  -  e 

v  jco  ) 


a  a+1  a 

j?+(xy-J±±  V  <  f0^r  > 


ro  +  i) 


A  -y®V+1 
1-e 


K  jco  J 


One-sided  power  functions: 


xa  x  >  0 

0,  otherwise 
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Symmetric  B-splines 


■  Symmetrization  in  Fourier  domain: 


a+ 1 

sin(W2) 

jo  ) 

<y/2 

0.5- 


-0.4 


#(*)-  F-1 


sin(<39/2) 
col  2 


a+ 1 


4 


3 


0 
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Properties 


Generic  notation  :  pa  for  either  /?“  (causal)  or  /?“  (symmetric) 

■  Equivalence  with  classical  B-splines 

P"(x)  with  a  =  n  (integers) 

P“(x)  with  a  =  2n  +  \  (odd  integers) 

Compact  support ! 

■  Decay 

-  ii  c 

Theorem  :  For  a  >  -1 ,  there  exists  a  constant  C  such  that  \/3a(x)\  <  . 

pc 

(U.  &  Blu,  SIAM  Rev,  2000) 

■  Convolution  property 

pax  *  pa2  _  pax+a2+ 1 

(/3a(-),j3a(--x))=tfa+l(x) 
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Riesz  basis 


j; /3a(x-k)}k  is  a  Riesz  basis  for  the  cardinal  fractional  splines 


Generic  B-spline  representation  of  a  fractional  spline 


s(x)  =  ^c[£]/?“(x  -  £) 

K  Zj 

Continuous-time  function 

Discrete  representation 

(fractional  spline) 

(digital  signal) 

Stable,  one-to-one  representation 

For  a  >-  j,  there  exist  two  constants  Aa>  0  and  Bn  <+oo  such  that 


Vce/2, 


< B  •  c 

«  II  II/, 
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Explicit  fractional  differentiation  formula 


■  Fractional  derivative  operators 

0  <  Foarier  >  ( jay 


Fractional  finite  difference  operator: 

£S  ^  Fourier.  ^  ^ 


Sketch  of  proof: 

0PaM)  < — > 


(j»1 


f  -\a+l 

\-eJ(0 


\  JO)  J 


'l 


JO)  J 
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Order  of  approximation 


Approximation  space  at  scale  a 

K  =  \sa O)  =  Xc(£)<pf“  ~k]:c(k)  e  l2 


k&Z 


Projection  operator 

V/'eZ2,  PJ  =  arg  min||/  -  sa  | 


V„ 


Sa  e  Va 


Order  of  approximation 
DEFINITION 


a  =  1 


1  2  3  4  5 


A  scaling  function  (p  has  order  of  approximation  y  iff 

V/sr,  \\f-P„j\<C-a'\f>\  =  0{a') 


B-splines  of  degree  a  have  order  of  approximation  y=a+1 
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Spline  reconstruction  of  a  CAT-scan 


Piecewise  constant 

Y  =1 


Cubic  spline 

y  =  4 
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Reproduction  of  polynomials 


More  fractals 


QuickTime™  and  a 
TIFF  (LZW)  decompressor 
are  needed  to  see  this  picture. 


Dali 


QuickTime™  and  a 
TIFF  (LZW)  decompressor 
are  needed  to  see  this  picture. 


Pollock 


QuickTime™  and  a 
TIFF  (Uncompressed)  decompressor 
are  needed  to  see  this  picture. 


Mandelbrot  meets  Mondrian 


FRACTIONAL  WAVELETS 


■  Basic  ingredients 

■  Constructing  fractional  wavelets 

■  Fractional  B-spline  wavelets 

■  Multi-scale  fractional  differentiation 

■  Adjustable  wavelet  properties 
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Scaling  function 


DEFINITION:  cp(x)  is  an  admissible  scaling  function  of  L2  iff: 

•  Riesz  basis  condition 


■  Two-scale  relation 


cp(x  /  2)  =  2^  h(k) cp(x  -  k ) 

i 

k  1  1 

keZ 

- ► 

■  Partition  of  unity 


2<p(*-fc)=  i 

keZ 
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From  scaling  functions  to  wavelets 


Wavelet  bases  of  Z_2  (Mallat-Meyer,  1989) 

For  any  given  admissible  scaling  function  of  4,  <p(x) ,  there  exits  a  wavelet 

Hx  /2)  =  2  YjS(k)(p(x  -  k ) 


k&Z 


such  that  the  family  of  functions 

[  j  r 


.V27 


-v 


x  —  2l  k 
2  j 


1  -1 


/eZ,kZ 

forms  of  Riesz  basis  of  4 . 


Constructive  approach:  perfect  reconstruction  filterbank 
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Constructing  fractional  wavelets 


Theorem  :  Let  (p{x)  be  the  L2 -stable  solution  (scaling  function)  of  the  two-scale  relation 
<p(x/2)  =  2^Jh(k)<p(x-k) 

ksZ 

Then  cp(x)  is  of  order  /  (fractional)  if  and  only  if 

v 

Q(z )  with  \Q(ejC0)\<cQ 

distributional  part 

spline  part 


Multi-scale  differentiator 

B-spline  factorization: 

Approximation  order: 

\J/(co)  oc(-y'(o)y,  to  — >  0 

<^> 

9=Pr‘*9o 

<^> 

ll/-^JL,  =0(a‘) 

11 


Vanishing  moments: 

J  xny/(x)dx  =  0,  n  =  0, •••[/-{] 


(Unser  &  Blu,  IEEE-SP ,  2003) 
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Binomial  refinement  filter 


■  Two-scale  relation 

/%(x/2)  =  2^h“(k)/%(x- k) 

k&Z 

■  Generalized  binomial  filter 


K(k)  =  - 


+1 


a  +  1 


Ha(z)  = 


1  +z 


2 


'ii)  rju  + 1) 

vvj  r(v  +  i)r(w- v  + 1) 


Example  of  linear  splines:  a=  1 


<2  +  1 
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Fractional  B-spline  wavelets 


JfceZ 


a  + 1 

V  n  ) 


pr\n  +  k-\)  pa+{x-k) 


g(k) 


QuickTime™  and  a 
Video  decompressor 
are  needed  to  see  this  picture. 


Remarkable  property 

Each  of  these  wavelets  generates  a  semi-orthogonal  Riesz  basis  of  L2 
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FFT-based  wavelet  algorithm 


■  Filterbank  algorithm 


(p(x / 2)  =  V2  y '',h(k)<p(x  -  k ) 

k&Z 


</Kx  12)  =  a/2  ^Yjg{k)(p{x  -  k) 

kzZ 


Click  for  demo 


(Blu  &  Unser,  ICASSP’ 2000) 
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Adjustable  wavelet  properties 


af 

/O) 

Wavelet  ( 

transform 

filtering 
detection 
feature  extraction 
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Transform  is  tunable  in  a  continuous  fashion  ! 


■  Order  of  differentiation:  y=a+1 
■  Whitening  of  fBMs,  fractals  .... 


■  Regularity 


■  Holder  continuity:  a 

■  Sobolev:  smax=  a+ 1/2 
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Wavelets  and  the  uncertainty  principle 


=  min 


|(x-x0M*)|L 


A(y=min 

(Do 


■  Heisenberg’s  uncertainty  relation 


Question:  are  there  such  wavelet  bases  ? 
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Localization  of  the  B-spline  wavelets 


Theorem 

The  B-spline  wavelets  converge  (in  Zp-norm)  to 
modulated  Gaussians  as  the  degree  goes  to  infinity  : 


lim  {/?»}=  C- 

a— »oo  v  J 


e~(x~xa  )2  /2cra 


\\m{y/+(x)}=  C' 

<2— » 00  V  J 


e~(x-x'a  filer 
Gaussian 


^  x  cos(<y0x  +  #a) 

sinusoid 


a 'a=B  cra  with  B  =  2.59 


Cubic  B-spline  wavelets: 
within  2%  of  the  uncertainty  limit ! 

(Unseret  a/.,  IEEE-IT,  1992) 


QuickTime™  and  a  Video  decompressor  are  needed  to  see  this  picture. 
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Are  there  wavelets 
in  my  brain  ? 


Fig.  2.  Similarity  between  the  receptive  field  of  simple  cortical 
cells  and  a  wavelet  basis  function,  (a)  Response  of  a  simple  X  cell 
from  a  monkey  visual  cortex  and  its  fined  Gabor  elementary  signal 
[26],  [67t  Fig.  (b)  Semi-orthogonal  cubic  E-spline  wavelet  and 
its  log-log  frequency  response  [100}. 


WAVELETS  IN  MEDICAL  IMAGING 


■  Survey  of  applications 

■  Analysis  of  functional  imaging  data  (fMRI) 
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Wavelets  in  medical  imaging: 
Survey  1991-1999 


References 

•  Unser  and  Aldroubi,  Proc  IEEE,  1996 

•  Laine,  Annual  Rev  Biomed  Eng,  2000 


Special  issue,  IEEE  Trans  Med  Im,  2003 


Image  pocessng  task  Application/ modify  Principal  Authors 


Image  cotpresson 

•  MR  Angdis94;  DeVcre  95 

•  Marrmogams  Manduca  95;  Wang  96; 

•  CT 

•  Angogams,  etc... 

Ritering 

Image  enhanamnt  Lajne94,95; 

•  Digital  radogams  Lu,  94;  Qan95; 

•MR  Guang97; 

•  Marrmogams  etc  . . . 

•  LungX-rays,  CT 

Dendsing  Weave-  9 1 ; 

•MR  Xu94;Coifman95j 

•  Ultrasound(speckle)  Abd^-IVfelek97;Laine98; 

•  SPECJ  Novak  98  99 

Feature  extraction 

Detection  ofrricio-calcficalions  Qian  95;  Yostida  94; 

•  Marrmogams  Strickland  96  Dhawan  96; 

Bacyu  96,  Heine  97V\fong98 

Texture  andysis anddassificetion  Barman  93  Lane94; Unser 

•  Ultrasound  95;  \Ne i  95  Yung95;  Bush 

•  CT  MR  97;  Mojsilovb97 

•  Mammograms 

Snakes  andactive  contours  ChuangKuo  96 

•  Ultrasound 

Wfavelet  encoding 

•  Magn&c  resonace  imagng  V\feava-l-baly  92; 

Parych  94,  96(3eman96; 
Shimzu96;Jian  97 

Image  reonstrudion 

•  Computer  tomogapty  Qson  93  94;  P^rn  94 

•  United  andp  data  V\folnut  93;  Dten^f  95 

•  Optcal  tomogapty  Sahner  96;Zhu97; 

.  SPECT  Kdaczyk  94;Raheja  99 

Statistical  data  analyss 

Fundiond  imagng  Ruttimann93,94,98; 

•  PET  Unser  95;  Feiher  99  Raz  99 

•fMR 

Multi -scale  Regstration 

Motion  correction  Unser  93;Thevena  95  95 

•  fMR,  angogapty  Kybic  99 

Multi -modality  imagng 

•  CT  PET,  MR 

3D  visualization 

•  CT  MR  Gross  95,  97|Vluraki  95; 

Kamath  98JHorbelt99 

QuickTime™  and  a  TIFF  (Uncompressed)  decompressor  are  needed  to  see  this  picture. 
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Functional  brain  imaging  by  fMRI 


BOLD  (Blood  Oxygenation  Level  Dependence) 

Basic  principle:  deoxygenated  blood  is  more  paramagnetic  than  oxygenated  blood 


Time  series 

A 
B 
A 
B 

A 

**  EPI  acquisition 

Matrix  size:  128  x  128  x  30  Pixels  x  68  measurements 
Resolution:  1.56  x  1.56  x  4  mm  x  6  seconds 


B:  Rest 
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Functional  brain  imaging  by  fMRI  (Cont’d) 


QuickTime™  and  a  decompressor  are  needed  to  see  this  picture. 


■  Main  problems: 

Small  signal  changes  (1-5%) 

Very  noisy  data  —  averaging 

■  Standard  solution 

■  Spatial  Gaussian  smoothing  (SPM) 
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On  the  fractal  nature  of  fMRI  data 


Brain:  courtesy  Jan  Kybic 


Log-Log  plot  of  spectral  density 


Fractal  dimension:  D  =  \  +  d  —  H  =  2.534  with  d  =  2  (topological  dimension) 
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Wavelet  analysis  of  fMRI 


■  Advantages  of  the  wavelet  transform  (Ruttiman  et  ai.,  ieee-tmi,  1998) 

Orthogonal  transformation  :  white  noise  -»  white  noise 

■  Decorrelates/whitens  fMRI  signal 

■  Data  compression 

■  Increased  signal-to-noise  ratio  (averaging  effect) 

Preserves  space  localization 


An  example:  auditory  stimulation 


Conclusion 


■  Fractional  splines 

■  Natural  extension  of  Schoenberg’s  polynomial  splines 

■  Stable,  convenient  B-spline  representation 

■  Most  polynomial  B-spline  properties  are  retained 

■  Intimate  link  with  fractional  calculus 

■  Elementary  building  blocks:  Green  functions  of  fractional  derivative  operators 

■  Efficient  digital-filter-based  solutions 

■  New  fractional  wavelets 

■  Multiresolution  bases  of  L2 
•  Fast  algorithm 

■  Tunable 

>  Regularity 

■  Localization 

■  Order  of  differentiation 

■  Optimal  for  the  processing  of  fractal-like  processes  (pre-whitening) 

■  Application  in  signal  and  image  processing 

■  Processing  of  fractal-like  signals 

■  Wavelet-based  processing  and  feature  extraction 
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Software  and  demos  at: 
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Extensions  (on-going  work) 


■  Richer  family:  alpha-tau  splines 


4  ( j4+\-j4~t 


[Blu  et  al.,  ICASSP’03] 


■  Multi-dimensional:  fractional  polyharmonic  splines 

■  Polyharmonic  smoothing  splines  [Tirosh  et  ai.,  icassp’04] 

■  Polyharmonic  wavelets  lVan  de  Vi|le  et  al  - under  review] 
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